{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "</center>\n",
    "Автор материала: программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Материал распространяется на условиях лицензии [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Тема 4. Линейные модели классификации и регрессии\n",
    "## <center>Часть 3. Наглядный пример регуляризации логистической регрессии"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "В 1 статье уже приводился пример того, как полиномиальные признаки позволяют линейным моделям строить нелинейные разделяющие поверхности. Покажем это в картинках.\n",
    "\n",
    "Посмотрим, как регуляризация влияет на качество классификации на наборе данных по тестированию микрочипов из курса Andrew Ng по машинному обучению. \n",
    "Будем использовать логистическую регрессию с полиномиальными признаками и варьировать параметр регуляризации C.\n",
    "Сначала посмотрим, как регуляризация влияет на разделяющую границу классификатора, интуитивно распознаем переобучение и недообучение.\n",
    "Потом численно установим близкий к оптимальному параметр регуляризации с помощью кросс-валидации (`cross-validation`) и  перебора по сетке (`GridSearch`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "\n",
    "# отключим всякие предупреждения Anaconda\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.linear_model import LogisticRegression, LogisticRegressionCV\n",
    "from sklearn.model_selection import GridSearchCV, StratifiedKFold, cross_val_score\n",
    "from sklearn.preprocessing import PolynomialFeatures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загружаем данные с помощью метода `read_csv` библиотеки `pandas`. В этом наборе данных для 118 микрочипов (объекты) указаны результаты двух тестов по контролю качества (два числовых признака) и сказано, пустили ли микрочип в производство. Признаки уже центрированы, то есть из всех значений вычтены средние по столбцам. Таким образом, \"среднему\" микрочипу соответствуют нулевые значения результатов тестов.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# загрузка данных\n",
    "data = pd.read_csv(\n",
    "    \"../../data/microchip_tests.txt\", header=None, names=(\"test1\", \"test2\", \"released\")\n",
    ")\n",
    "# информация о наборе данных\n",
    "data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на первые и последние 5 строк."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.tail(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сохраним обучающую выборку и метки целевого класса в отдельных массивах NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = data.iloc[:, :2].values\n",
    "y = data.iloc[:, 2].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Отобразим данные. Красный цвет соответствует бракованным чипам, зеленый – нормальным.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X[y == 1, 0], X[y == 1, 1], c=\"green\", label=\"Выпущен\")\n",
    "plt.scatter(X[y == 0, 0], X[y == 0, 1], c=\"red\", label=\"Бракован\")\n",
    "plt.xlabel(\"Тест 1\")\n",
    "plt.ylabel(\"Тест 2\")\n",
    "plt.title(\"2 теста микрочипов\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Определяем функцию для отображения разделяющей кривой классификатора"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_boundary(clf, X, y, grid_step=0.01, poly_featurizer=None):\n",
    "    x_min, x_max = X[:, 0].min() - 0.1, X[:, 0].max() + 0.1\n",
    "    y_min, y_max = X[:, 1].min() - 0.1, X[:, 1].max() + 0.1\n",
    "    xx, yy = np.meshgrid(\n",
    "        np.arange(x_min, x_max, grid_step), np.arange(y_min, y_max, grid_step)\n",
    "    )\n",
    "\n",
    "    # каждой точке в сетке [x_min, m_max]x[y_min, y_max]\n",
    "    # ставим в соответствие свой цвет\n",
    "    Z = clf.predict(poly_featurizer.transform(np.c_[xx.ravel(), yy.ravel()]))\n",
    "    Z = Z.reshape(xx.shape)\n",
    "    plt.contour(xx, yy, Z, cmap=plt.cm.Paired)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Полиномиальными признаками до степени $d$ для двух переменных $x_1$ и $x_2$ мы называем следующие:\n",
    "\n",
    "$$\\large \\{x_1^d, x_1^{d-1}x_2, \\ldots x_2^d\\} =  \\{x_1^ix_2^j\\}_{i+j=d, i,j \\in \\mathbb{N}}$$\n",
    "\n",
    "Например, для $d=3$ это будут следующие признаки:\n",
    "\n",
    "$$\\large 1, x_1, x_2,  x_1^2, x_1x_2, x_2^2, x_1^3, x_1^2x_2, x_1x_2^2, x_2^3$$\n",
    "\n",
    "Нарисовав треугольник Пифагора, Вы сообразите, сколько таких признаков будет для $d=4,5...$ и вообще для любого $d$.\n",
    "Попросту говоря, таких признаков экспоненциально много, и строить, скажем, для 100 признаков полиномиальные степени 10 может оказаться затратно (а более того, и не нужно). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Создадим объект `sklearn`, который добавит в матрицу $X$ полиномиальные признаки вплоть до степени 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = PolynomialFeatures(degree=7)\n",
    "X_poly = poly.fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_poly.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обучим логистическую регрессию с параметром регуляризации $C = 10^{-2}$. Изобразим разделяющую границу.\n",
    "Также проверим долю правильных ответов классификатора на обучающей выборке. Видим, что регуляризация оказалась \n",
    "слишком сильной, и модель \"недообучилась\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = 1e-2\n",
    "logit = LogisticRegression(C=C, n_jobs=-1, random_state=17)\n",
    "logit.fit(X_poly, y)\n",
    "\n",
    "plot_boundary(logit, X, y, grid_step=0.01, poly_featurizer=poly)\n",
    "\n",
    "plt.scatter(X[y == 1, 0], X[y == 1, 1], c=\"green\", label=\"Выпущен\")\n",
    "plt.scatter(X[y == 0, 0], X[y == 0, 1], c=\"red\", label=\"Бракован\")\n",
    "plt.xlabel(\"Тест 1\")\n",
    "plt.ylabel(\"Тест 2\")\n",
    "plt.title(\"2 теста микрочипов. Логит с C=0.01\")\n",
    "plt.legend()\n",
    "\n",
    "print(\n",
    "    \"Доля правильных ответов классификатора на обучающей выборке:\",\n",
    "    round(logit.score(X_poly, y), 3),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Увеличим $C$ до 1. Тем самым мы *ослабляем* регуляризацию, теперь в решении значния весов логистической регрессии могут оказаться больше (по модулю), чем в прошлом случае. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = 1\n",
    "logit = LogisticRegression(C=C, n_jobs=-1, random_state=17)\n",
    "logit.fit(X_poly, y)\n",
    "\n",
    "plot_boundary(logit, X, y, grid_step=0.005, poly_featurizer=poly)\n",
    "\n",
    "plt.scatter(X[y == 1, 0], X[y == 1, 1], c=\"green\", label=\"Выпущен\")\n",
    "plt.scatter(X[y == 0, 0], X[y == 0, 1], c=\"red\", label=\"Бракован\")\n",
    "plt.xlabel(\"Тест 1\")\n",
    "plt.ylabel(\"Тест 2\")\n",
    "plt.title(\"2 теста микрочипов. Логит с C=1\")\n",
    "plt.legend()\n",
    "\n",
    "print(\n",
    "    \"Доля правильных ответов классификатора на обучающей выборке:\",\n",
    "    round(logit.score(X_poly, y), 3),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Еще увеличим $C$ – до 10 тысяч. Теперь регуляризации явно недостаточно, и мы наблюдаем переобучение. Можно заметить, что в прошлом случае (при $C$=1 и \"гладкой\" границе) доля правильных ответов модели на обучающей выборке не намного ниже, чем в 3 случае, зато на новой выборке, можно себе представить, 2 модель сработает намного лучше. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = 1e4\n",
    "logit = LogisticRegression(C=C, n_jobs=-1, random_state=17)\n",
    "logit.fit(X_poly, y)\n",
    "\n",
    "plot_boundary(logit, X, y, grid_step=0.005, poly_featurizer=poly)\n",
    "\n",
    "plt.scatter(X[y == 1, 0], X[y == 1, 1], c=\"green\", label=\"Выпущен\")\n",
    "plt.scatter(X[y == 0, 0], X[y == 0, 1], c=\"red\", label=\"Бракован\")\n",
    "plt.xlabel(\"Тест 1\")\n",
    "plt.ylabel(\"Тест 2\")\n",
    "plt.title(\"2 теста микрочипов. Логит с C=10k\")\n",
    "plt.legend()\n",
    "\n",
    "print(\n",
    "    \"Доля правильных ответов классификатора на обучающей выборке:\",\n",
    "    round(logit.score(X_poly, y), 3),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Чтоб обсудить результаты, перепишем формулу для функционала, который оптимизируется в логистической регрессии, в таком виде:\n",
    "$$J(X,y,w) = \\mathcal{L} + \\frac{1}{C}||w||^2,$$\n",
    "\n",
    "где\n",
    " - $\\mathcal{L}$ – логистическая функция потерь, просуммированная по всей выборке\n",
    " - $C$ – обратный коэффициент регуляризации (тот самый $C$ в `sklearn`-реализации `LogisticRegression`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Промежуточные выводы**:\n",
    " - чем больше параметр $C$, тем более сложные зависимости в данных может восстанавливать модель (интуитивно $C$ соответствует \"сложности\" модели (model capacity))\n",
    " - если регуляризация слишком сильная (малые значения $C$), то решением задачи минимизации логистической функции потерь может оказаться то, когда многие веса занулились или стали слишком малыми. Еще говорят, что модель недостаточно \"штрафуется\" за ошибки (то есть в функционале $J$ \"перевешивает\" сумма квадратов весов, а ошибка $\\mathcal{L}$ может быть относительно большой). В таком случае модель окажется *недообученной* (1 случай)\n",
    " - наоборот, если регуляризация слишком слабая (большие значения $C$), то решением задачи оптимизации может стать вектор $w$ с большими по модулю  компонентами. В таком случае больший вклад в оптимизируемый функционал $J$ имеет  $\\mathcal{L}$ и, вольно выражаясь, модель слишком \"боится\" ошибиться на объектах обучающей выборки, поэтому окажется *переобученной* (3 случай)\n",
    " - то, какое значение $C$ выбрать, сама логистическая регрессия \"не поймет\" (или еще говорят \"не выучит\"), то есть это не может быть определено решением оптимизационной задачи, которой является логистическая регрессия (в отличие от весов $w$). Так же точно, дерево решений не может \"само понять\", какое ограничение на глубину выбрать (за один процесс обучения). Поэтому $C$ – это *гиперпараметр* модели, который настраивается на кросс-валидации, как и *max_depth* для дерева."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Настройка параметра регуляризации**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь найдем оптимальное (в данном примере) значение параметра регуляризации $C$. Сделать это можно с помощью `LogisticRegressionCV` – перебора параметров по сетке с последующей кросс-валидацией. Этот класс создан специально для логистической регрессии (для нее известны эффективные алгоритмы перебора параметров), для произвольной модели мы бы использовали `GridSearchCV`, `RandomizedSearchCV` или, например, специальные алгоритмы оптимизации гиперпараметров, реализованные в `hyperopt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=17)\n",
    "\n",
    "c_values = np.logspace(-2, 3, 500)\n",
    "\n",
    "logit_searcher = LogisticRegressionCV(Cs=c_values, cv=skf, verbose=1, n_jobs=-1)\n",
    "logit_searcher.fit(X_poly, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logit_searcher.C_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, как качество модели (доля правильных ответов на обучающей и валидационной выборках) меняется при изменении гиперпараметра $C$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(c_values, np.mean(logit_searcher.scores_[1], axis=0))\n",
    "plt.xlabel(\"C\")\n",
    "plt.ylabel(\"Mean CV-accuracy\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выделим участок с \"лучшими\" значениями C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(c_values, np.mean(logit_searcher.scores_[1], axis=0))\n",
    "plt.xlabel(\"C\")\n",
    "plt.ylabel(\"Mean CV-accuracy\")\n",
    "plt.xlim((0, 10));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Такие кривые называются *валидационными*, и в `sklearn` для них их построения есть специальные методы."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
